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Superconducting pairing of electrons in nanoscale metallic particles with discrete energy levels 
and a fixed number of electrons is described by the reduced BCS model Hamiltonian. We show 
that this model is integrable by the algebraic Bethe ansatz. The eigenstates, spectrum, conserved 
operators, integrals of motion, and norms of wave functions are obtained. Furthermore, the quantum 
inverse problem is solved, meaning that form factors and correlation functions can be explicitly 
evaluated. Closed form expressions are given for the form factors and correlation functions that 
describe superconducting pairing. 
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Due to recent advances in nanotechnology it has be- 
come possible to fabricate and characterise individual 
metallic grains with dimensions as small as a few nanome- 
ters They are sufhciently small that the spacing, d, of 
the discrete energy levels can be determined. A particu- 
larly interesting question concerns whether superconduc- 
tivity can occur in a grain with d comparable to A, the 
energy gap in a bulk system. If c? <C A, the supercon- 
ducting correlations are well-described by a mean-field 
solution to the reduced pairing Hamiltonian (equation 
(|l|) below) due to Bardeen, Cooper, and Schrieffer (BCS) 
in the grand canonical ensemble with a variable number 
of electrons. However, if c? ~ A recent numerical calcu- 
lations have shown that when the number of electrons is 
fixed (as in the canonical ensemble) the superconducting 
fluctuations become large and approximate treatments 
become unreliable [0j2|. Thus, exact calculations of phys- 
ical quantities are highly desirable. It has only recently 
been appreciated that the exact eigenstates and spectrum 
of the BCS model were found in the 1960's by Richard- 
son, in the context of nuclear physics The model 
has subsequently been found to have a rich mathemati- 
cal structure: it is integrable (i.e., has a complete set of 
conserved operators) Q|, has a connection to conformal 
field theory , and is related to Gaudin's inhomogeneous 
spin- 1/2 models [|-§. 

In this Letter we show how the BCS model can be 
solved using the algebraic Bethe ansatz (ABA) method. 
This result can be deduced from the observation that the 
conserved operators obtained in 0| were also obtained in 
1^ via the ABA, but in another context. However, the ap- 
proach we adopt here is slightly different from ||^ , which 
facilitates the solution of the quantum inverse problem 
|pO|-p^ to explicitly evaluate form factors (i.e., one point 
functions) and correlation functions. This completes the 
agenda recently set out by Amico, Falci, and Fazio 
We also readily obtain known results for eigenstates, the 
spectrum, and conserved operators. Our treatment is 
also applicable to superconductivity in fermionic atom 



traps |l|jT| and can also be extended to a solvable model 
for condensate fragmentation in boson systems [l^ . 

The Hamiltonian for the reduced BCS model consists 
of a kinetic energy term and an interaction term which 
describes the attraction between electrons in time re- 
versed states, 
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where j = 1, • • • , labels a shell of doubly degenerate 
single particle energy levels with energies and Cj^ the 
annihilation operators; a = -\-,— labels the degenerate 
time reversed states; g denotes the BCS pairing couphng 
constant. Using the pseudo-spin realization of electron 



the BCS Hamiltonian (Q) becomes (up to 
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The R matrix. An essential ingredient of the ABA, 
which follows from the Quantum Inverse Scattering 
Method (QISM), is the construction of the R matrix solv- 
ing the quantum Yang-Baxter equation, 

i?12(wi - U2)^13("l - '"3)i?23(w2 " ^g) 

= -R23(W2 - W3)i?i3(ui - U3)Ri2{ui - U2) 

where the Uj are spectral parameters. Here Rjk denotes 
the matrix oiyV ®V ®V (where V is the two-dimensional 
Hilbert space on which the pseudo-spin operators act) 
acting on the j-th and /c-th spaces and as an identity 
on the remaining space. The R matrix may be viewed 
as the structural constants for the Yang-Baxter algebra 
generated by the monodromy matrix T(u), 

Rl2{ui - U2)Ti{ui)T2{u2) = T2[u2)Ti{ui)Ri2{ui - U2). 

(3) 
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There are two kinds of realisations of the Yang-Baxter 
algebra which are relevant to our construction. One is 
operator- valued given by the R matrix Roj{u) and the 
other is a c-number representation G which does not de- 
pend on the spectral parameter u. In the latter case, we 
have [R{u), G® G] = 0. The comultiplication behind the 
Yang-Baxter algebra allows us to construct a representa- 
tion of the monodromy matrix through 

Tq{u) = GoRoniu ~ en) ■ ■ •Goi?oi('" - ei)- 

Defining the transfer matrix via t(u) = troToiu) it fol- 
lows that [t{u),t{v)] — for all values of the parameters 
u, V. If the R matrix possesses the regularity property 
Rjk{0) — Pjk with P being the permutation operator, 
then it is easily verified that 



^(cj) = GjRjj-i{ej - dj-i) ■ ■ ■ GjRji{e 
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xGjRjn{ej ~ £12) • • ■ Gji?j.,+i (e^ - ej+i)Gj. 

Let us now assume that the R matrix is quasi-classical, 
i.e., it admits a series expansion R(u) = / + rir[u) + • • •, 
for an appropriate parameter rj. If we can also choose G 
such that G — 1 + rjT + ■ ■ ■, then the expansion of t{ej) 
in terms of 77 takes the form. 
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An immediate consequence from the commutativity of 
the transfer matrics is [Tj,Tfe] = 0. Therefore an inte- 
grable model is obtained by taking the set {ri} as the 
conserved operators and a Hamiltonian given as a func- 
tion of the Tj. 

We apply the procedure described above to the su{2) 
invariant R matrix R{u) = h(u)I + c(u)P, with en- 
tries that are rational functions: b{u) = u/(u -\- rj) and 
c(u) = ri/(u + rj). Note that the regularity property 
i?(0) = P is present. For this case, we can choose G 
as any element of the su{2) algebra. We claim that the 
BCS model corresponds to the special choice 



Gj = e^pi~2riS^/gn). 
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This can be viewed as a generalized inhomogeneous six- 
vertex model. Expanding this and the R matrix to first 
order in 77 and substituting in (||) we find from that 
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where we have discarded a constant term. These oper- 
ators are the isotropic Gaudin Hamiltonians in a non- 
uniform magnetic field Their relevance to the spin 
realization of the BCS model (Q) is that the latter is ex- 
pressible (up to a constant) as 
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Although the above expressions for the consereved opera- 
tors only applies to the case when all tj 's are distinct, our 
construction can be adapted to accommodate the cases 
when some of e^'s are the same. 

Algebraic Bethe ansatz. In the ABA, the integrals of 
motion are obtained by finding the eigenfunctions of the 
transfer matrix which is given by the trace of the mon- 
odromy matrix. The monodromy matrix is written in the 
form 



T{u) 



A{u) B{u) 
G{u) D{u) 



which is the quantum equivalent of the scattering coef- 
ficients of the classical inverse scattering problem. Then 
from the Yang-Baxter algebra (^, we may derive the 
fundamental commutation relations (FCR) between the 
entries of the monodromy matrix. Choosing the state 
|0) — (E'jLiI T)i as the pseudovacuum, then we have 
the pseudovacuum eigenvalues a{u) and d{u) of A{u) and 
D{u): a{u) = eyip{-'q / g) , d{u) = eyip{ril g)Y{^b{u - tj). 
Following the standard procedure pO| , p^ , we choose the 
Bethe state 
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Then we may derive the off-shell Bethe ansatz equations 
using the FCR following which, in the quasiclas- 

sical limit, takes the form 
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where 
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In (|^) we defined i/'a by 



J = l 



Imposing fa ~ 0, one immediately sees that -0 becomes 
the eigenvector of the conserved operator Tj with Xj as 
the eigenvalue. The constraint = is then equivalent 
to Richardson's equations 0] , 
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Here L ~ N may be interpreted as the number of time- 
reversed pairs of electrons. The energy eigenvalue of the 
Hamiltonian (0) is 
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Scalar products and norms. Directly evaluating the 
norms of Bethe wave functions can be tedious, if not 
impossible. However, using the QISM they can be repre- 
sented as determinants . Since this representation 
only depends on the R matrix, the derivation presented 
previously for different models can be readily applied to 
our (generalized) inhomogeneous six-vertex model. In 
the QISM construction, the determinant representation 
for scalar products 

N N 

{0\X{C{wp)X{B{v^)\Q) 

/3=1 a=l 

play a crucial role; especially, when one of the sets of 
parameters, for example {va}, is a solution of the Bethe 
equations pO| , pl] , pl| . In the quasiclassical limit, the lead- 
ing term of the scalar product for the inhomogeneous 
six- vertex model gives rise to the scalar product 



with K = HjLi Gj = exp{-2riJ2j=i S^/gQ). The above 
construction is one of our main results. The appearance 
of the powers of K arises from the c-number matrix re- 
alisation of the Yang-Baxter algebra G which is peculiar 
to our construction. Following one can obtain the 
representation of the correlation functions in terms of 
pseudovacuum eigenvalues a{u) and d{u). 

Form factors. For the BCS model the pair correlator 



is of particular interest ||l|,^ . (We use the notation that 
(x) = {vir ■ -vnWW, - ■ ■ ,vn)/{vi, - ■ -vnW, - ■ ■ ,vn) for 
any operator x)- can be interpreted as the prob- 

ability enhancement of finding a pair of electrons in 
level TO, instead of two uncorrelated electrons. (It is 
zero for g = 0). In the pseudo-spin representation 
Cl = {S-^S+){S+S^) = 1I^-{S^„,)\ In general, form 
factors such as 
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where the matrix elements of J are given by 



can be calculated for the generalized inhomogeneous six- 
^20)vertex model. In the quasiclassical limit, they reduce to 
the form factors of the BCS model. 
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Here are a solution to Richardson's equations (j^), 

whereas \wp} are arbitrary parameters. Richardson's ex- 
pression 122] for the square of the norm of the Bethe state 
follows from ([l0| ) and (^l]) by taking the limit Wa ^ Va. 

Solution of the quantum inverse problem. In order to 
calculate the form factors and correlation functions, we 
need to solve this problem for the generalized inhomo- 
geneous six-vertex model. This then allows the recon- 
struction of local quantum spin operators in terms of the 
quantum monodromy matrix. A general procedure for 
doing this has recently been presented |jll|,|l^ for the so- 
called fundamental models. In our case, where the model 
is not fundamental, we find 
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Above, T is the N x N matrix obtained from T by delet- 
ing the last row and column and replacing iV + 1 by A'^ 
in the matrix elements. Here we assume that both {wq} 
and {wh} are solutions to Richardson's Bethe equations 
(^. However, the results are still valid for if only 
{wb} satisfy the Bethe equations. 

Correlation functions. We find that the correlation 
functions of the BCS model take the same form as the un- 
derlying su{2) spin 1/2 Gaudin model, with the param- 
eters Vj satisfying Richardson's Bethe ansatz equations 
instead of Gaudin's ones. Here we present explicitly 
the two-point correlation function 
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Further applications. Our work is also relevant to 
proposals to observe BCS superconductivity in gases of 
fermionic atoms such as spin-polarised ^Li |lj] . Quantum 
degeneracy of ^Li at temperatures of about 240 nK has 
recently been observed in an atom trap with frequencies, 
Lu ^ I kHz ||l^ , corresponding to an energy level spacing 
of the order of 10~^^ eV. The estimated BCS transition 
temperature is of the order of 20 nK , corresponding 
to an energy gap of the order of 4 x 10^^^ cV. Hence, 
these systems are in a regime where the physics consid- 
ered here will be important. 

Dukelsky and Schuck |l^ recently introduced a solv- 
able model for condensate fragmentation in finite boson 
systems. The model they solved follows from the con- 
struction used above when the Yang-Baxter algebra is 
realized in terms of the generators of the Lie algebra 
szt(l, 1). The model also provides a new mechanism for 
the enhancement of sd dominance in interacting boson 
models in the context of nuclear physics p6| . 

This work was supported by the Australian Research 
Council. 



Here the hat denotes that the corresponding parameter 
is not present in the set. Since {wa} is a solution of 
the Bethe equations, {wi, • • • , wn\S^\vi, • • • , Wq,, • • • , vn) 
is the form factor given before, while 

{wi, - ■ ■ ,wn\s~s:;;^\vi, - ■ ■ ,VN-2) = 

U0=iiwf3 - em)(wff - £«) dctjvT(m, n, {wp}, {«„}) 
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In (|l|) TO ^ n is assumed, with the convention that it 
is zero when m = n. The above results constitute the 
building blocks of the Penrose-Onsager-Yang off-diagonal 
longe-range order (ODLRO) parameter Aqd Q, 



^OD = ^ ^(5+S'„). 



(15) 



The small grain behavior of this parameter and its con- 
nection with the pair correlator (12) was recently dis- 
cussed in 1251. 
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